function [Y] = MenoSRMain(TSTART, DELT, NSTEP, FREQ, HMEAS)

global NTHE NPHIM icalc
global GX GY GZ
global SPINM AXAs AYAs AZAs espin spinr
global T1mx T1my T1mz T2mx T2my T2mz
global JEXG RR AR BR ZFS temp rwid nrval rr2 fr2 fnon
global GXB GYB GZB AXBs AYBs AZBs
global blim T1Npar T1Nper
global relwid gxh gyh gzh T1h T2h frh

% Main function to calculate baseline saturation curve

if nrval == 1
  fr1 = 1.0;
else
  fr1 = 1.0 - fr2;
end
if fnon > 0.0
  fr1 = fr1 - fnon;
end
if frh > 0.0
  fr1 = fr1 - frh;
end
if fr1 < 0.0
  h = errordlg(['Fraction of conformation 1 is negative.' ...
    '  Please re-enter values.'], 'ERROR');
  uiwait(h);
  return;
end

% TODO: What the heck is this?
if icalc == 2
  nometal(y,dat);
  % change parameters
end

% value of beta is in cm-1/G, w0 is in cm-1
beta = 0.46686E-4;
w0 = 0.033356*FREQ;
nlinr = 2*spinr + 1;  % number of radical hyperfine lines
% set up limits for orientation selection
blow = HMEAS - 0.5*blim;
bhi = HMEAS + 0.5*blim;
% limdist is number of points to define distribution of relaxation rates
limdist = 9;
if rwid < 0.01
  limdist = 1;
end
% limr is number of points to define distribution of relaxation rates
limr = 9;
if relwid < 0.01
  limr = 1;
end
% zero array to hold sr curve
Y(NSTEP) = 0.0;
probt = 0.0;
TEND = TSTART + (NSTEP - 1)*DELT;
% CONVERT ANGLES TO RADIANS
ARR = AR*0.01745;
BRR = BR*0.01745;
% CALCULATE AVERAGE G FOR METAL
GAVM = 0.333*(GX+GY+GZ);
% CONVERT J IN GAUSS TO J IN RAD/SEC
JEX = single(JEXG*(GAVM + 2.0059)*4.398E6);
% disp(['GAVM: ' num2str(GAVM)])
% disp(['JEXG: ' num2str(JEXG)])
% disp(['JEX: ' num2str(JEX)])
% CALCULATE ORIENTATION-INDEPENDENT PART OF DIPOLE INTERACTION
% 2.04e10 = 1.0e24 * (beta squared)/(hbar *4) where 1.0e24 is conversion
% of r from Angstroms to cm and 4 in denominator is to match factor
% in the "B" term
% calculation of "part" uses rr for conformation 1
PART = 2.0059*2.04E10/(RR^3);
% CALCULATE resonant frequency for radical IN RAD/SEC, this is
% for center of spectrum, not including hyperfine here.
WNIT = FREQ*2.0E9*3.1416;
% parameters for metal nuclear spin
NLINE = 2*SPINM + 1;
SP2 = SPINM*(SPINM+1);
% NOW CALCULATE FACTORS THAT ARE INDEPENDENT OF
% ORIENTATION
GX2 = GX*GX;
GY2 = GY*GY;
GZ2 = GZ*GZ;
GXB2 = GXB*GXB;
GYB2 = GYB*GYB;
GZB2 = GZB*GZB;
gxh2 = gxh*gxh;
gyh2 = gyh*gyh;
gzh2 = gzh*gzh;
% calculate terms needed for orientation dependence of T1N
% disp(T1Nper)
% disp(T1Npar)
del = 1/T1Nper - 1/T1Npar;
if spinr ~= 0.0
  % CONVERT radical HYPERFINE to CM-1
  AXB = 1.0E-4*AXBs;
  AYB = 1.0E-4*AYBs;
  AZB = 1.0E-4*AZBs;
  AXB2 = AXB*AXB;
  AYB2 = AYB*AYB;
  AZB2 = AZB*AZB;
  ANONB2 = (AYB2-AXB2)*(AYB2-AXB2);
  AXYB2 = AXB2*AYB2;
end
if SPINM ~= 0.0
  % CONVERT metal HYPERFINE FROM 10-4 TO RAD/SEC
  AXA = 1.885E7*AXAs;
  AYA = 1.885E7*AYAs;
  AZA = 1.885E7*AZAs;
  AXA2 = AXA*AXA;
  AYA2 = AYA*AYA;
  AZA2 = AZA*AZA;
  ANON2 = (AYA2-AXA2)*(AYA2-AXA2);
  AXYA2 = AXA2*AYA2;
end
GXY2 = GX2*GY2;
GXZ2 = GX2*GZ2;
GYZ2 = GY2*GZ2;
if espin > 0.5
  % calculate 2D/kT including 1.986e-16 for conversion of cm-1 to erg
  % and k = 1.38e-16 erg/K. Alternatively think of this as dividing
  % ZFS by 0.695 to convert to Kelvin
  sep3 = 2.878*ZFS/temp;
  pop3 = exp(-1.0*sep3);
  if espin == 1.5
    P1= 1.0/(1.0+pop3);
    P3= pop3/(1.0+pop3);
    W2D=1.885E11*2*ZFS;
  end
  if espin == 2.5
    sep5 = 3.0*sep3;
    pop5 = exp(-1.0*sep5);
    P1 = 1.0/(1.0+pop3+pop5);
    P3 = pop3/(1.0+pop3+pop5);
    P5 = pop5/(1.0+pop3+pop5);
    % to convert ZFS to rad/sec, multiply by 2*pi*speed of light
    W2D = 1.885E11*2*ZFS;
    W4D = 2*W2D;
    %     disp(['wnit' num2str(WNIT)]);
    %     disp(['w2d' num2str(W2D)]);
    %     disp(['w4d' num2str(W4D)]);
  end
  %   disp(['p1' num2str(p1)]);
  %   disp(['p3' num2str(p3)]);
  %   disp(['p5' num2str(p5)]);
end
% set up parameters for powder averaging
% HALF OF THE STEPS ARE IN EQUAL INCREMENTS OF THETA UP
% TO 43 DEGREES AND THEN REST ARE IN EQUAL STEPS OF
% COS THETA
NTHEH = floor(NTHE/2);
% CALCULATE INCREMENT IN THETA UP TO 43 DEGREES
STEP1 = 0.7505/NTHEH;
% CALCULATE INCREMENT IN COS THETA FOR 43 - 90 DEGREES
STEP2 = 0.7314/NTHEH;
% START OF THETA LOOP
THEP = 0.7505;
hWbar = waitbar(0, 'Calculating Saturation Recovery Curve');
for IT = 1:NTHE
%   disp(['THETA CNTR: ' num2str(IT)]);
  waitbar(IT/NTHE, hWbar);
  if IT <= NTHEH
    THE = (2*IT-1)*STEP1/2;
    COSTHE = cos(THE);
    DELTH = STEP1;
  else
    ITT = IT - NTHEH;
    COSTHE = 0.7314 - (2*ITT-1)*STEP2/2;
    THE = acos(COSTHE);
    DELTH = THE-THEP;
    if ITT == 1
      DELTH = 2*(THE-0.7505);
    end
    THEP = THE;
  end
  CSTHE2 = COSTHE*COSTHE;
  SNTHE2 = 1-CSTHE2;
  SINTHE = sqrt(SNTHE2);
  ST2CT2 = SNTHE2*CSTHE2;
  % calculate R1NO for this orientation
  R1NO = 1/T1Npar + SINTHE*SNTHE2*del;
  % CALCULATE NUMBER OF PHI INCREMENTS FOR THIS VALUE OF THETA
  NPHI = floor(NPHIM*SINTHE+1);
  STEP = 6.2832/NPHI;
  % NEED TO SUM OVER 4 OCTANTS, DO THIS WITH THETA FROM 0 TO 90 AND
  % PHI FROM 0 TO 360
  % START OF PHI LOOP
  for IP = 1:NPHI
    % disp(['PHI counter: ' num2str(IP)]);
    PHI = (2*IP-1)*STEP/2;
    if NPHI == 1
      PHI = 0.0;
    end
    COSPHI = cos(PHI);
    SINPHI = sin(PHI);
    CSPHI2 = COSPHI*COSPHI;
    SNPHI2 = SINPHI*SINPHI;
    % NOW EVALUATE TERMS INVOLVING THETA AND PHI for metal
    ST2CP2 = SNTHE2*CSPHI2;
    ST2SP2 = SNTHE2*SNPHI2;
    FGX = GX2*ST2CP2;
    FGY = GY2*ST2SP2;
    FGZ = GZ2*CSTHE2;
    GSQ = FGX + FGY + FGZ;
    G = sqrt(GSQ);
    % WMETC IS FREQUENCY of metal resonance, before considering nuclear
    % hyperfine splitting
    % 8.799E6 = 2*pi*beta/h
    WMETC = G*HMEAS*8.799E6;
    % analogous terms for hemichrome
    if frh > 0.0
      gsqh = gxh2*ST2CP2 + gyh2*ST2SP2 + gzh2*CSTHE2;
      gh = sqrt(gsqh);
      wmeth = gh*HMEAS*8.799E6;
    else
      gsqh = 0;
      gh = 0;
      wmeth = 0;
    end
    % WPRO accounts only for solid angle
    WPRO = SINTHE*DELTH*STEP;
    WPROB = WPRO/(NLINE*nlinr);
    GPERP2 = GX2*CSPHI2+GY2*SNPHI2;
    % calculate metal relaxation times at this orientation,
    % assume orientation dependence can be done analogous to A values
    % for distributions, these are central values
    R1MS = ((1/T1mx)^2*FGX+(1/T1my)^2*FGY+(1/T1mz)^2*FGZ)/GSQ;
    T1MC = sqrt(1/R1MS);
    R2MS = ((1/T2mx)^2*FGX+(1/T2my)^2*FGY+(1/T2mz)^2*FGZ)/GSQ;
    T2GC = sqrt(1/R2MS);
    %     if IP == 1
    %     disp(['T1mc' num2str(T1MC)]);
    %     disp(['T2gc' num2str(T2GC)]);
    %     end
    if SPINM == 0
      AK = 0.0;
      R4R5 = 0.0;
      R1R2R3 = 0.0;
    else
      AK2 = (AXA2*FGX + AYA2*FGY + AZA2*FGZ)/GSQ;
      AK = sqrt(AK2);
      APERP2 = (AXA2*GX2*CSPHI2+AYA2*GY2*SNPHI2)/GPERP2;
      % CALCULATE 2ND ORDER TERMS FOR METAL NUCLEUS
      R12 = AZA2*APERP2/AK2;
      R22 = AXYA2/APERP2;
      % CALCULATE PARTS OF R32 AND R52 WHICH ARE THE SAME
      RP2 =(ANON2*GX2*GY2*SNPHI2*CSPHI2)/(AK2*GSQ*GPERP2);
      R32 = (RP2*AZA2*GZ2*CSTHE2)/(APERP2*GPERP2);
      ADIF = APERP2-AZA2;
      R42=(ST2CT2*GZ2*GPERP2*ADIF/(AK2*GSQ*GSQ))*ADIF;
      R52 = RP2*SNTHE2;
      R1R2R3 = (R12+R22+R32)/(4*WMETC);
      R4R5 = (R42+R52)/(2*WMETC);
    end
    % now do analogous terms for radical
    FGXB = GXB2*ST2CP2;
    FGYB = GYB2*ST2SP2;
    FGZB = GZB2*CSTHE2;
    GSQB = FGXB + FGYB + FGZB;
    GB = sqrt(GSQB);
    BETAG = beta*GB;
    BETAG2 = BETAG*BETAG;
    H0 = w0/BETAG;
    GPERPB2 = GXB2*CSPHI2+GYB2*SNPHI2;
    if spinr == 0.0
      AKB = 0.0;
    else
      AKB2 = (AXB2*FGXB + AYB2*FGYB + AZB2*FGZB)/GSQB;
      AKB = sqrt(AKB2);
      APERPB2 = (AXB2*GXB2*CSPHI2+AYB2*GYB2*SNPHI2)/GPERPB2;
      % CALCULATE 2ND ORDER TERMS FOR radical NUCLEUS
      R12B = AZB2*APERPB2/AKB2;
      R22B = AXYB2/APERPB2;
      % CALCULATE PARTS OF R32 AND R52 WHICH ARE THE SAME
      RP2B =(ANONB2*GXB2*GYB2*SNPHI2*CSPHI2)/(AKB2*GSQB*GPERPB2);
      R32B = (RP2B*AZB2*GZB2*CSTHE2)/(APERPB2*GPERPB2);
      ADIFB = APERPB2-AZB2;
      R42B=(ST2CT2*GZB2*GPERPB2*ADIFB/(AKB2*GSQB*GSQB))*ADIFB;
      R52B = RP2B*SNTHE2;
      R1R2R3B = (R12B+R22B+R32B)/(4*BETAG2);
      R4R5B = (R42B+R52B)/(2*BETAG2);
    end
    % COSDELR IS COSINE OF ANGLE BETWEEN MAGNETIC FIELD AND INTERSPIN VECTOR
    COSDELR = sin(ARR)*SINTHE*cos(BRR-PHI)+cos(ARR)*COSTHE;
    SINDELR2 = 1.0 - COSDELR*COSDELR;
    % these are terms that depend on r^-6 because part depends on r^-3
    % Can simply scale csq and esq for different r but can't do that for
    % bsq so move bsq inside of loop over r's
    CSQ = GSQ*((6.0*PART*COSDELR)^2)*SINDELR2;
    ESQ = GSQ*((3.0*PART*SINDELR2)^2);
    % SAME, BUT FOR HEMICHROME
    CSQH = gsqh*((6.0*PART*COSDELR)^2)*SINDELR2;
    ESQH = gsqh*((3.0*PART*SINDELR2)^2);
    %     if IP == 1
    %     disp(['cosdelr' num2str(COSDELR)]);
    % %     disp(['bsq' num2str(BSQ)]);
    %     disp(['csq' num2str(CSQ)]);
    %     disp(['esq' num2str(ESQ)]);
    %     end
    % START OF LOOP OVER RADICAL SPIN STATES
    for ispn = 1:nlinr
      spinn = spinr-ispn+1;
      spn2 = spinn*spinn;
      HEST = H0-AKB*spinn/BETAG;
      HRES = HEST-spn2*R4R5B/HEST - R1R2R3B*(spinr*(spinr+1)-spn2)/HEST;
      if (HRES > blow) && (HRES < bhi)
        % add in contribution to SR curve if orientation selected
        % loop over conformations
        for ir = 1:nrval
          if ir == 1
            wr = fr1;
            rmid = RR;
          else
            wr = fr2;
            rmid = rr2;
          end
          % loop over r distribution
          for ik = 1:limdist
            rcal = rmid-(5-ik)*rwid/2.0;
            % rra and rrat and corrections to part and part^2
            rra = (RR/rcal)^3;
            rrat = (RR/rcal)^6;
            wtd = wr*WPROB*exp(-0.173*(5-ik)^2);
            BSQ = (-0.5*JEX-G*PART*rra*(1.0-3.0*COSDELR*COSDELR))^2;
            BSQH = (-0.5*JEX-gh*PART*rra*(1.0-3.0*COSDELR*COSDELR))^2;
            % start of loop over r distribution
            for mr = 1:limr
              if limr == 1
                T1m = T1MC;
                T2g = T2GC;
                wtda = wtd;
              else
                % calculation done in steps of 0.25*relwid,
                % which is width at half
                % height of distribution
                if mr < 5
                  T1m = T1MC*(1+((5-mr)/4.0)*relwid);
                  T2g = T2GC*(1+((5-mr)/4.0)*relwid);
                else
                  T1m = T1MC/(1+((mr-5)/4.0)*relwid);
                  T2g = T2GC/(1+((mr-5)/4.0)*relwid);
                end
                wtda = wtd/4.231*exp(-0.173*(5-mr)^2);
              end
              if espin == 0.5
                INT2 = (2.0*BSQ);
                DELC = CSQ*rrat*T1m/(1.0+(WNIT*T1m)^2);
                DELE = 2.0*ESQ*rrat*T2g/(1.0+((WMETC+WNIT)*T2g)^2);
              end
              if espin == 1.5
                DELB = BSQ*(6*(P3+P1)*T2g/(1.0+((W2D-WNIT)*T2g)^2)...
                  +8*P1*T2g/(1.0+((WMETC-WNIT)*T2g)^2+BSQ*T1m*T2g));
                DELC = CSQ*rrat*((9*P3+P1)*T1m/(1+(WNIT*T1m)^2)+...
                  3*(P3+P1)*T1m/(1+(W2D*T1m)^2)+4*P1*T1m/(1+(WMETC*T1m)^2));
                DELE = ESQ*rrat*(6*(P3+P1)*T2g/(1.0+((W2D+WNIT)*T2g)^2)+...
                  8*P1*T2g/(1.0+((WMETC+WNIT)*T2g)^2));
              end
              if espin == 2.5
                DELB = BSQ*(10*(P5+P3)*T2g/(1.0+((W4D-WNIT)*T2g)^2)+...
                  16*(P3+P1)*T2g/(1.0+((W2D-WNIT)*T2g)^2)+...
                  18*P1*T2g/(1.0+((WMETC-WNIT)*T2g)^2+BSQ*T1m*T2g));
                DELC = CSQ*rrat*((25*P5+9*P3+4*P1)*T1m/(1+(WNIT*T1m)^2)+...
                  5*(P5+P3)*T1m/(1+(W4D*T1m)^2)+8*(P3+P1)*T1m/...
                  (1+(W2D*T1m)^2)+9*P1*T1m/(1+(WMETC*T1m)^2));
                DELE = ESQ*rrat*(10*(P5+P3)*T2g/(1+((W4D+WNIT)*T2g)^2)+...
                  16*(P3+P1)*T2g/(1.0+((W2D+WNIT)*T2g)^2)+18*P1*T2g/...
                  (1.0+((WMETC+WNIT)*T2g)^2));
              end
              if frh > 0.0
                INT2 = (2.0*BSQH);
                DELOM = wmeth - WNIT;
                DELBh = (INT2*T2h)/(1.0+(DELOM*T2h)^2+INT2*T1h*T2h);
                DELCh = CSQH*rrat*T1h/(1.0+(WNIT*T1h)^2);
                DELEh = 2.0*ESQH*rrat*T2h/(1.0+((wmeth+WNIT)*T2h)^2);
                WR1NOH = R1NO+(DELBh+DELCh+DELEh);
              end
              % Start of Loop Over Metal Spin States
              for ISP = 1:NLINE
                probt = probt + wtd;
                if espin == 0.5
                  % assume we don't have nuclear hyperfine
                  % for S = 5/2 (Fe) or S = 3/2 (CR)
                  SPIN = SPINM-ISP+1;
                  SPIN2 = SPIN*SPIN;
                  WMET = WMETC+AK*SPIN+SPIN2*R4R5+R1R2R3*(SP2-SPIN2);
                  DELOM = WMET - WNIT;
                  DELB = (INT2*T2g)/(1.0+(DELOM*T2g)^2+INT2*T1m*T2g);
                end
                WR1NO = R1NO+(DELB+DELC+DELE);
                %  if ISP == 1
                %    disp(['R1NO, WR1NO' num2str(R1NO) ' ' num2str(WR1NO)]);
                %    disp(['DELOM, DELB' num2str(DELOM) ' ' num2str(DELB)]);
                %    disp(['DELB, DELC, DELE ' num2str(DELB) ' ' num2str(DELC) ' ' num2str(DELE)]);
                %  end
                for IND = 1:NSTEP
                  TIME = TSTART +(IND-1)*DELT;
                  Y(IND) = Y(IND) - wtda*exp(-1.0*TIME*WR1NO);
                end
              end % end of loop over metal spin states
            end % loop over distribution in relaxation times
            % assume that hemichromes do not have distribution
            % of relaxation times
          end % end of loop over distance distribution
        end % end of loop over conformations
        % addition of contribution from non-interacting radical needs to be
        % inside the if-block for orientation selection
        if fnon > 0.0
          for IND = 1:NSTEP
            TIME = TSTART + (IND-1)*DELT;
            Y(IND) = Y(IND) - (WPRO/nlinr)*fnon*exp(-1.0*TIME*R1NO);
          end
          % wpro accounts for solid angle, nlinr is number of radical lines
          % wprob and wtd account for solid angle, number of metal hyperfine
          % lines, and number of radical lines
        end
        if frh > 0.0
          for IND = 1:NSTEP
            TIME = TSTART + (IND-1)*DELT;
            Y(IND) = Y(IND) - (WPRO/nlinr)*frh*exp(-1.0*TIME*WR1NOH);
          end
        end
        % the following end is for the end of orientation selection
      end
    end % end of loop over radical spin states
  end % end of PHI loop
end % end of theta loop
close(hWbar);
% disp(['Total Probt ' num2str(probt)]);
% disp(['Y: ' num2str(Y)]);
% disp(['nstep, tstart, delt' num2str(NSTEP) ' ' num2str(TSTART) ' ' num2str(DELT)]);
if probt ~= 0
  % Convert time to micro_seconds
  TSTARTM = TSTART*1E6;
  DELTM = DELT*1E6;
%   disp(['first, last ' num2str(Y(1)) ' ' num2str(Y(NSTEP))])
%   disp(['TSTAT, DELTM : ' num2str(TSTARTM) ' ' num2str(DELTM)])
  % eprplt(Y, dat, NSTEP, LEGEND, NLINES, TSTARTM, DELTM, 2);
  % setting final parameter in call as 2, identifies this as menosr plot
else
  h = errordlg('No contribution at this field', 'Error');
  uiwait(h);
end
return
end